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I. INTRODUCTION 
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In this paper, we examine in detail the principal branches of solutions that arise in vector dis- 
crete models with nonlinear inter-component coupling and four wave mixing. The relevant four 
branches of solutions consist of two single mode branches (transverse electric and transverse mag- 
netic) and two mixed mode branches, involving both components (linearly polarized and elliptically 
polarized). These solutions are obtained explicitly and their stability is analyzed completely in the 
anti-continuum limit (where the nodes of the lattice are uncoupled), illustrating the supercritical 
pitchfork nature of the bifurcations that give rise to the latter two, respectively, from the former 
two. Then the branches are continued for finite coupling constructing a full two-parameter numer- 
ical bifurcation diagram of their existence. Relevant stability ranges and instability regimes are 
highlighted and, whenever unstable, the solutions are dynamically evolved through direct compu- 
tations to monitor the development of the corresponding instabilities. Direct connections to the 
Oh earlier experimental work of Meier et al. [Phys. Rev. Lett. 91, 143907 (2003)] that motivated the 

present work are given. 
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^"J ■ Recently, nonlinear Hamiltonian lattice dynamical systems with a large number of degrees of freedom have become 
a focal point for a variety of application areas Several diverse physical contexts in which such models (typically 
discrete in space and continuous in the evolution variable) arise are (i) spatial dynamics of optical beams in coupled 
waveguide arrays arising in nonlinear optics 0] , (ii) temporal evolution of Bose-Einstein condensates (BECs) in optical 
lattices in soft-condensed matter physics and (iii) DNA double strand in biophysics Q . 
y—i • A prototypical model, applicable to different degrees of approximation, in all of the above contexts is the so-called 
| discrete nonlinear Schrodinger (DNLS) equation. This model was first proposed in [j| and implemented for the first 
time experimentally in Refs. 0,0- A systematic presentation of the experimental results in optical systems is given 
in a recent work || ; for BEC-related experiments, see 0, . In this model, self-localized excitations (discrete soiitons) 
are possible as a result of the interplay between the Kerr nonlinearity and discrete linear coupling. Many properties of 
optical discrete spatial soiitons have been systematically explored in theory and experiment, including generalizations 
ly-j , to diffraction management [H1IDJ> diffraction-managed soiitons an d soliton transport and gating |l3l ITiL HHfill . 

On the other hand, a topic that has received considerably less attention has been the study of vector analogs of 
""^3 \ the DNLS equation. A number of studies have addressed the existence and stability of diverse families of solitary 
waves/localized states in the DNLS equation (see 0,^3,01)- These issues have also been studied for both cubic and 
quadratic nonlinearities [Tsl |20| in one and two dimensions |l7l l2ll |2^ . However, the first experimental realization 
for such a system occurred rather recently (see also the very recent paper |24[). This work sheds new light 
on the theoretical investigations of this topic, in that the relevant model for the corresponding AlGaAs waveguide 
array experiment included four wave mixing (FWM) terms that were rarely included in previous theoretical studies. 
Subsequent studies (see H^,|2^|) addressed a number of specific properties of the model system put forward in [23| . 
Such properties related to switching, instability-induced amplification, modulational instability as well as an analysis 
of the energy barrier between stable discrete solitary waves centered on a lattice site versus two lattice sites. 

The focus of the present study is to extend the analysis of the vector NLS lattice system incorporating FWM which 
is relevant to the experiments of p3 | . More specifically, we provide a systematic analysis of the existence and stability 
requirements of the principal solitary wave modes of the system. In so doing, we have found that the concept of the 
anti-continuum limit (in which the coupling between waveguides is absent) provides a powerful tool for the study of the 
existence and stability concerning these solutions. Subsequently, continuations in the coupling parameter are used to 
construct a full two-parameter bifurcation diagram that gives the complete picture of stability for both one-component 
as well as "mixed" mode solutions (involving both components). We note that all of the results obtained via our 
dimensionless model can be directly placed in the context of the experimental results obtained in Refs. [23ll2^ |. 
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The remainder of the paper is organized in the following way. In section II, we present the theoretical model of 
interest as well as proposing a setup that allows for analysis of the existence and stability of the model's solutions. In 
section III, we fully analyze the relevant solutions and their stability in the anti-continuum limit (i.e., zero coupling 
between adjacent waveguides). In section IV, we proceed numerically to determine solutions to the model for nonzero 
coupling and observe how the bifurcation points and stability features are modified for this case. The instability 
regions obtained in section IV are then monitored through direct numerical experiments in section V. Finally, in 
section VI, we summarize our findings and present our conclusions. 



II. MODEL AND SETUP 

Following |2^], we will use the dimensionless model 

ia n = -a n - e(a n+ i + a n _i) - (|a„| 2 + A\b n \ 2 ) a n - Bb 2 n a* n (1) 

ib n = b n -e(b n+1 +b n ^ 1 )- (\b n \ 2 + A\a n \ 2 ) b n - Ba 2 n b* n . (2) 

In these equations, a n and b n are the appropriately normalized, slowly-varying, complex field envelopes for the 
transverse electric (TE) and transverse magnetic (TM) polarized waves respectively. The constants A and B are 
respectively associated with the cross-phase modulation (XPM) and four-wave mixing (FWM) and were evaluated 
in [23 to be approximately equal to A ~ 1 and B ~ 1/2. Wherever possible, the results obtained in this analysis 
will be given for general A and B, even though the above values have been used in the specifics of our numerical 
computations. The overdot in Eqs. CJ-J2J denotes differentiation with respect to the evolution variable (which is z 
in this case) and the * denotes the complex conjugate. For the properties of the waveguide array and incident light 
used in the experiment (see |2^| for details), the dimensionless power 

p = J2(K\ 2 + K\ 2 ) (3) 

n 

is connected with its dimensional analog Pd (measured in Watts) according to Pd ~ 56. 4P. Furthermore, the dimen- 
sionless coupling e used in [23j was e ~ 0.921. This is crucial since it serves to highlight differences between the results 
obtained in the computational analysis of |23j and the generalization of the corresponding results presented herein. 
More generally, we take e to be a free parameter in our analysis. We note that e is related to the dimensional coupling 
constant, k, given in |23| . via e = 2.741k. Here, k is measured in mm -1 . 

We look for stationary solutions in the form: a n = a n e lqz and b n = b n e tqz and subsequently dropping the tildes, we 
obtain the corresponding stationary equations: 

(q - l)o„ - e K+i + o n _i) - (K| 2 + A\b n \ 2 ) a n - Bb 2 n a* n = (4) 

(q + 1)6„ - e (b n+1 + bn-!) - (\b n \ 2 + A|a„| 2 ) b n - Ba 2 n b* n = 0. (5) 

In this context, q is an additional free parameter, namely the propagation constant (q is also used as a free parameter 
for the computations presented in |23j). Our parameter q is related to the corresponding dimensional propagation 
constant, S, given in |23j . by S — 0.3648<7. Here, S is measured in mm -1 . 

It is worth pointing out that for Fig. 1 in Ref. (2^|, the range of parameters used in their computations appear to 
correspond to the case of q > 3. This is verified by our analysis below. Again, this is highlighted to account for the 
disparities of some of the conclusions obtained in |23j | with the ones their generalizations presented in the following 
sections. 

To perform linear stability computations, provided that an exact solution (a°,6°) of the stationary Eqs. (0J-© 
has been obtained (either analytically in the anti-continuum limit, or numerically for e ^ 0), we use the perturbation 
ansatz: 

a n = a° n + S(c n e~ l " z + d n e l ^ z ^j (6) 
b n = b° l + 5(f n e-^ z +g n e^ z ) . (7) 

In this ansatz, S is a (small) formal expansion parameter and u> is the eigenfrequency of the perturbation eigenvector 
(c„, d* , f n , <7*) T (where the superscript denotes the transpose). Then, upon lengthy but straightforward calculation, 
one obtains the linearized matrix eigenvalue problem to 0(<5). These dynamical equations have the form: 
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The N x N (where N is the size of the lattice) blocks of the linearization matrix are given by: 



L n = ( q -l)-e(A 2 + 2)-2K\ 2 -A\b°f (8) 

L12 = -Kf-B{b° n ? (9) 

L13 = -Aa° n (b n y-2B(a n rb n (10) 

L 14 = -Aa° n b° n (11) 

£21 = -L* a (12) 

L22 = -L n (13) 

L23 = ~L* U (14) 

L 2 4 = —L\ 3 (15) 

L 31 = L* 13 (16) 

£32 = £14 (17) 

£33 = ((7+l)- £ (A 2 + 2)-2|6 r J 2 -AK| 2 (18) 

£34 = -(6°) 2 -S«) 2 (19) 

L 41 = -L* u (20) 

£42 = -Lis (21) 

L43 = ^Lg 4 (22) 

L44 = — L33 (23) 



In the above, the short-hand notation (A 2 + 2)z n = 2 n +i + z n-ii invoking the concept of the discrete Laplacian A 2 . 
With this linearized matrix eigenvalue problem, we can establish a criterion for the existence and stability of solutions 
to this system by examination of our system at the anti-continuum limit (i.e., when e = 0). 



III. THE ANTI-CONTINUUM LIMIT 



A. Existence 



We begin by examining eqs. (@J — JSJ) for the case e = 0. Setting a n = r n e l6n and b n — s n e l ^ n in this case, one 
finds: 

(q - 1) - [rl + Asl) Bs 2 n e 2i ^- 9 ^ = (24) 

(q + 1) - (si + Arl) Br*e- M <*»-«»> = 0. (25) 

From the above, we infer that for solutions to exist in this limit it is necessary for 2(9 n — <j> n ) to be an integer multiple 
of 7T, hence, we obtain that 

$n-K = k\ (26) 

with k G Z. 

The simplest possible solutions are the ones that involve only one of the two branches and were hence termed TE 
and TM modes respectively in |2^|. The TE solution of Eqs. I|24|) - (|25|l has the form (in the present limit) 



r n =±Vg-l s„=0 (27) 
and exists only for q > 1. On the other hand, the TM mode features 

r„ = o s „ = ±VgTT, (28) 
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and is only present for q > — 1. 

Also, there are two possible mixed mode solutions permitted by the quantization condition of cq. I|2()|) . allowing 
the exponential of Eqs. (|24|l - 1)25(1 to be +1 or —1 respectively. The first one (e 2 ^ 9 ™ - ^™) = 1) was characterized as a 
linearly polarized (LP) branch in [23|, involving in-phase contributions from both the TE and TM components. In 
this case, the linear system of Eqs. <|24ll - H25[) has the general solution: 



r i_, (A + B)( q + l)- iq -l) 

n ~ V (A + B) 2 - 1 [ J) 



(A + B)(q-l)-(q + l) 
Sn = ± V (A + B) 2 -l ' (30) 

If (A + B) 2 > 1 (which is the case in the experimentally relevant setting), this branch only exists for (A + B)(q + 1) > 
(q — 1) and (A + B)(q— 1) > {q+l) [the sign of the two above inequalities needs to be reversed for existence conditions 
in the case of (A + B) 2 < 1]. Among the two conditions, in the present setting, the second one is the most "stringent" 
for the case A = 2B = 1, which yields the constraint q > 5 (while the first condition requires for the same parameters 
?>-5). 

The second possibility for a mixed mode (stemming from e 2 ^ e " ^ n > = —1 in Eqs. I|24|l - H25|l ') yields the, so-called, 
elliptically polarized mode (EP) [2j|, involving a it/ 2 phase shift between the TE and TM components. The relevant 
amplitudes in this case are: 



(A-B)(q + l)-(q-l) 
rn = ±sl {A _ B y_ l (3D 



(A-B)(q~l)-(q + l) 
Sn = ± \/ {A-B) 2 -l ' (32) 

The condition for this branch to exist if (A - B) 2 < 1 is q — 1 > (A — B)(q + 1) and q+l > (A — B){q — 1) (once 
again the signs should be reversed if (A — B) 2 > 1). In this case, the first condition is more constraining than the 
second, imposing for our special case of interest herein (A = 2B = 1) q > 3 (while the second condition only requires 
<Z>-3). 



B. Stability 

We can now turn to the study of the stability of the corresponding branches via the anti-continuum limit (e = 0). The 
major advantage of this limit is that all the blocks of the linearization matrix L now become diagonal. Furthermore, 
the structures that we consider, per the branches presented above, involve the excitation of a single site (all other sites 
are inert) with amplitudes corresponding to the appropriate TE, TM, LP or EP modes. The interesting byproduct of 
this formulation is that all inert sites end up yielding purely diagonal entries in the stability matrix equal to (q — 1) 
in the first (2N-1) x (2N-1) blocks and (q + 1) in the second set of (2N-l)x (2N-1) blocks. Hence the matrix will have 
2N — 1 eigenfrequencies equal to (q — 1) and equally many eigenfrcquencies of (q + 1). The remaining 4x4 matrix 
will correspond to the single excited site of the lattice and its eigenvalues will yield the eigenfrequencies pertaining to 
the corresponding branch of solutions. Hence, by solving this considerably simpler eigenvalue problem, we can infer 
the full stability picture for e =/= from the anti-continuum limit. This, we believe, is a major advantage in using this 
approach. 

We note that due to the fact that the power is conserved for our system (see previous section) , two of the relevant 
eigenvalues of this 4x4 problem should be identically 0. This is confirmed in the computations presented here. 

Now, taking the case for the "TE" mode in eqs.© — (|23|l . we determine the remaining pair of eigenvalues in this 
case: 

UJTE = ±y/(q + I) 2 + (A 2 - B 2 )(q - l) 2 - 2A(q 2 - 1). (33) 

Examining our model for the experimental case, A = 1, B = 1/2 and q > 0, we find that this eigenfrequency is real 
for q < 5, while it is imaginary for q > 5 (hence, implying the presence of an instability, since the corresponding 
eigenvalue A = iui becomes real and exponential growth occurs oc e Az ). Indeed, we can show that for the TE mode, 
stability is ensured for 1 < q < 5. 



FIG. 1: The full bifurcation diagram for e = 0, illustrating the two supercritical pitchfork bifurcations at q = 3 and q = 5. 
These give rise to the formation of the EP branch (the former) and the LP branch (the latter). The blue lines illustrate stable 
branches, while the red lines show the unstable ones. 

For the TM mode, the corresponding expression can be obtained as: 

ojtm = ±V(q - l) 2 + (A 2 - B 2 )(q + l) 2 - 2A(qi - 1), (34) 

and for the case of interest, stability is ensured for — 1 < q < 3, while instability ensues (due to a real eigenvalue) for 
q > 3. 

For the LP mode, the expression for the corresponding eigenfrequency is too cumbersome to provide in analytic 
form (even though such a form has been obtained) for general A and B. Hence, we will restrict ourselves to the case 
of A = 2B = 1, where the relevant eigenfrequency reads: 

uj LP = i^y^TT^ (35) 
5 

implying stability for q > 5 (which coincides with the branch's existence region). 
Finally, for the EP mode, the expression similarly reads: 

uj EP ±^y/gt~9, (36) 

implying, in principle, stability for q > 3 (again, coinciding with the solutions domain of existence); see, however, the 
discussion below for this branch. 

From the above observations, one can put together the complete bifurcation picture of the four branches and their 
corresponding stability in the anti-continuum limit (again, we present this here for concreteness for A — 2B = 1). 
More specifically, the TE branch exists for q > 1 and is stable for 1 < q < 5. For q > 5, the branch is destabilized 
as a new branch emerges, namely the LP branch, through a pitchfork bifurcation; notice that the TM component of 
this branch, per Eq. 1)30(1 is exactly zero at q = 5, hence it directly bifurcates from the TE branch. The two branches 
of this super-critical pitchfork correspond to the two signs of s n in Eq. (|3(J[I . The bifurcating branch "inherits" the 
stability of the TE branch for all larger values of q, while the latter branch remains unstable thereafter. 

In a similar development, the TM branch exists for q > —1 and is stable in the interval —1 < q < 3. However, at 
q = 3, a new branch (in fact, a pair of branches) bifurcates acquiring non-zero r„, beyond the bifurcation point, as per 
Eq. I|31|) . This is accompanied by the destabilization of the TM branch (due to a real eigenvalue) and the apparent 
stability of the ensuing EP branch for all values of q > 3. 

Fig. ^ shows how the two supercritical pitchfork bifurcations are generated as a function of the coordinates 
(r n ,s n ,q). We can now turn on the coupling and see how the existence and stability of each branch is affected by 
varying q and e. 
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IV. FINITE COUPLING 

One of the key features that will help us in our analysis of the finite coupling case consists of determining where 
the bands of the continuous spectrum of the linearized problem occurs. This is obtained analytically using the plane 
wave ansatz a n ~ ^(kn-uz) and ^ _ e i(kn-u>z)_ p\ug ging our ansatz into the coupled linearized problem (essentially, 
using eqs. Q — © and ignoring the nonlinear terms), the resulting dispersion relations are 

ui = (q-1) -2ecos(fc) (37) 
u) = {q + 1) - 2ecos(fc). (38) 

Hence the continuous spectrum (in the infinite lattice case) will consist of the frequency intervals [q — 1 — 2e, q — 1 + 2e] 
and [q + 1 — 2e, q + 1 + 2e] . As may be expected, for e = 0, these intervals degenerate to two single points (ui = q — 1 and 
uj = q + 1), in consonance with the discussion of the previous section. But perhaps more importantly, for q — 1 < 2e 
(equivalently for e > (q — l)/2), the continuous spectrum branch will be crossing the origin, leading to the collision 
of the eigenvalues with their mirror symmetric opposites, that will in turn lead to instabilities. Hence, we need only 
consider couplings in the interval e € [0, (q — l)/2]. 

We now numerically examine each of the main four branches, using one-parameter continuation in the above- 
mentioned interval of e's for different values of q and thus constructing a two-parameter bifurcation diagram for each 
branch. Notice that the solutions are numerically obtained for each pair of (q, e) by means of a fixed point iteration. 
Once convergent to within a pre-set tolerance (10~ 8 typically for the computations presented herein), the procedure 
is followed by numerical linear stability analysis, obtaining the eigenfrequencies for the linearization matrix L. 

A. TE branch 

The continuation of the TE branch is detailed in Fig. It turns out that especially in the case of this branch, 
solutions cannot be obtained for e > (q— 1)/2, i.e., the branch terminates at that point with its amplitude going to zero 
at this critical point. Within its region of existence, the branch has a domain of stability and one of instability. The 
point of separation between the two in the anti-continuum limit, studied previously, was the critical point of q = 5. 
For e ^ 0, the separatrix curve is shown in Fig. [2] and can be well approximated by the curve € TE ~ (4v / 2/5)-y / 9 _ -~ 5. 
Hence for q < 5, the solution is stable for all values of e in its range of existence (0 < e < (q — l)/2), while for q > 5, 
the solution is only stable for e^ E < e < (q — l)/2 and unstable (due to a real eigenvalue pair) for < e < e^E- 

B. TM branch 

The TM branch is considerably more complicated, structurally, than its TE counterpart. Firstly, it does not 
disappear beyond the critical e — (q — l)/2, however, it does become unstable as predicted previously, hence we will, 
once again, restrict ourselves to this parameter range. Furthermore, similarly to the TE branch case, there is an d TM 
below which the branch is always unstable, whereas for e > e^ Ml the branch may be stable. At the anti-continuum 
limit, the critical point for the instability is q = 3, as discussed previously; for q > 3, the critical point is obtained 
numerically in Fig. |3] It can be well approximated (close to q = 3) by e^ M ~ (9/10)y/q — 3. 

However, within the range of potential stability (0 < e < (q — l)/2 for q < 3, and e TM < e < (q — l)/2 for q > 3), 
we observe an additional large region of instability in the two parameter bifurcation diagram of Fig. |3| due to a 
complex quartet of eigenvalues. This instability appears to "emanate" from the point with (q, e) = (2.2, 0) in the 
anti-continuum limit, and to linearly expand its range as e increases. Hence, we will try to understand it at the level 
of e = 0. What happens at the value of q = 2.2 in the e = limit is that the point spectrum eigenfrequency of Eq. H34f) 
"collides" with the continuous spectrum point of concentration , co rresponding to u> = q— 1. However, the eigenvector 
of this eigenvalue [34( has a Krein signature (see e.g. [27], |2^, |2!j for relevant definitions of this signature) opposite 
to that of the continuous band at oj = q — 1. Hence, according to Krein theory |30(, the resulting collision leads to 
the formation of a quartet of eigenvalues emerging in the complex plane and, in turn, implying the instability of the 
TM configuration. As e is varied from 0, the continuous spectrum band grows linearly in e, hence the corresponding 
interval of g's, where this instability is present (due to the collision with the opposite Krein sign eigenvalue) also grows 
at the same rate. Along the same vein, it is worth pointing out that the line of this instability threshold and that of 
e = (q— l)/2 are parallel. 
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FIG. 2: TE branch: The top left set of panels shows the details of the TE branch as a function of e for q — 3. The topmost 
panel shows the norm of the branch as a function of e, while the left and right panels show the solution profile and its spectral 
plane for e = 0.1 and e = 0.9 respectively. The top right set of panels shows similar details, but for q — 5.5. The top right panel 
in the right case, which is unstable for small e, shows the most unstable eigenvalue real part. The middle and bottom panels 
show the (unstable) solution and its stability for e = 0.5 and the (stable) one for e = 1.5 respectively. The bottom panel shows 
the two parameter bifurcation diagram of the coupling e as a function of q (for details see text). All the relevant existence and 
stability regimes have been accordingly labeled. The top left and top right set of panels can be considered as two "vertical 
cuts" across this bifurcation diagram for q — 3 and 5.5 respectively. 

C. LP branch 

The mixed LP branch, involving in-phase contributions of the TE and TM has been continued for various values 
of q > 5, as a function of e. Recall that this branch emerges through a supercritical pitchfork as q is varied for for 
fixed e (see Fig. In an exactly analogous way, if e is varied for a fixed q > 5, the branch appears to terminate at 
e = e^E through a subcritical pitchfork. Hence, the LP branch only arises for < e < e^E an d <Z > 5 and it is stable 
throughout its interval of existence. For e > e^ E , this branch "degenerates" back into the stable (for such e's) TE 
mode. This is clearly illustrated in Fig. 01 that shows the continuation of the LP branch as a function of e for a fixed 
q = 5.5. The top panels show the norms of the two components illustrating that the second component is absent for 
e > e^ E — 0.81 in this case. The middle and bottom left panels show a case for e = 0.5 < e^ E where a genuine (and 
stable) TE solution exists, while the corresponding right panels are for e = 1 > ej, E) P as t the critical point where the 
LP solution degenerates into a TE mode. 

D. EP branch 

The case of the EP branch, analyzed in Fig. is somewhat analogous to that of the LP branch. For fixed q close 
to (and larger than) 3 and varying e, the branch exists and is stable for < e < e^ M . For larger values of e, the EP 
mode degenerates into a TM mode, as is shown in the top left panel of Fig. 0for q = 3.5 (cf. with the analogous for 
the LP case of Fig. 0}. 
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FIG. 3: TM branch: the figures are similar as in Fig. El but now for the TM branch. The top left set of panels shows 
the continuation in e for q = 2.5, along with the special cases of e = 0.1 (stable) and e = 0.5 (unstable due to quartet of 
eigenvalues) in the left and right middle and bottom panels. The top right set of panels shows the q = 3.5 (unstable) case with 
the middle/bottom panels corresponding to e = 0.5 (unstable due to real eigenvalue pair)/e = 0.65 (stable) respectively. 



However, for q > 3.62, this phenomenology appears to change and an expanding (for increasing q) interval of 
oscillatory instability within the range of existence of the EP branch (0 < e < e^ M ) seems to arise. To rationalize this 
feature, we again turn to the anti-continuum limit. Just as in the case of the (TM) branch from which it originates, 
the EP branch has a point spectrum eigenfrequency given by Eq. 1|36|) which has a negative Krein signature and upon 
collision with the continuous spectrum band of eigenfrequencies will give rise to an instability. Setting the frequency 
of Eq. I|36|) equal to q — 1, we obtain that this collision occurs at q = 9. For lower values of q, this "collision" will 
occur for a finite (nonzero) interval of values of e -see e.g. the top right panels of Fig. [SJfor q — 4.5-, which, in turn, 
yields the region of oscillatory instability of the EP mode, clearly separated from its region of stability by a dashed 
line in the bottom panel of Fig. |SJ 

In conclusion, the EP solutions exist in the regime where the TM mode is unstable (cf. bottom panels of Figs. 
and[SJ), illustrated by the curve originating from the point (g, e) = (3,0) in the anti-continuum limit. However, the 
EP solution, as is shown in the top panels of the figure, may (top left) or may not (top right) necessarily be stable in 
its region of existence due to the appearance of an oscillatory instability for sufficiently large q. 



V. DYNAMICAL EVOLUTION OF UNSTABLE BRANCHES 



We now examine a series of different instability scenarios pertaining to the branches of solutions discussed above. In 
so doing, we solve Eqs. Q-l© using a fourth-order Runge-Kutta algorithm. To ensure the validity of our numerical 
results, we have monitored the conservation of the power of Eq. and of the Hamiltonian (energy) of the system: 

H = - («an+i + a n a* n+1 + b* n b n+1 + b n b* n+1 ) + \a n \ 2 - \b n \ 2 + i(|a„| 4 + \b n \ 4 + A\a n \ 2 \b n \ 2 ) + ^(al(b* n ) 2 + b 2 n (a*) 2 ))(39 
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FIG. 4: LP branch: The top left and top right panel shows the dependence on e of the respective powers of the two components 
Pi and Pi (P — Pi + P2), for q — 5.5. The middle and bottom panel panel show the profile and the linear stability of the LP 
mode while it is still present (e = 0.5; left panels) and of the TE mode that it degenerates to, beyond e^E ( e = 1; right panels). 



The power is conserved to numerical precision, while the Hamiltonian is conserved to 1 part in 10 8 . In all cases, 150 
spatial sites are used in our computations. 

Our first example, shown in Fig. |HJ illustrates the case of q = 5.5 and e = 0.5 for the TE mode. This is a case 
representative of the region of parameter space (cf. Fig. |2J), where the TE mode is unstable due to a real eigenvalue 
pair (while the LP mode is stable). The solution actually coincides with the steady TE branch until approximately 
time 10. The configuration then deviates from the TE branch and subsequently develops a breathing oscillation 
around a stable LP state. 

The next three cases are concerned with the different scenarios for the instability of the TM mode. Fig. [7| with 
q = 3.5 and e = 0.5, is representative of the region of parameter space where the TM mode is unstable due to a real 
eigenvalue pair, while the EP mode is stable. In this case, As in the previous case, the solution coincides with the 
steady TM mode up until approximately time 10. It is then observed that for short times the configuration deviates 
from the TM branch and oscillates around the (stable) EP configuration. However, for longer times, the system 
appears to wander away from the latter configuration and into a breathing state, where the two components have 
roughly equal power and exchange a small fraction of it (almost) periodically. 

The second instability scenario of the TM branch, shown in Fig. |SJ involves the case of an eigenvalue quartet as 
e.g. for q — 2.5 and e = 0.5 (cf. Fig. In this case, we observe the configuration departing from the unstable TM 
steady state and resulting into a breathing solution involving both components. In this and the remaining cases, the 
departure from the steady state occurs later than in the first two cases. 

Next, we examine the evolution of the instability for the TM mode in the presence, also, of continuous spectrum 
instabilities for q = 2.5 and e = 1 (cf. Fig. |3J). In this case a different phenomenology occurs. The mode is completely 
destroyed in favor of small amplitude excitations. 

Finally, we examine the evolution of the instability for an elliptically polarized mixed mode case within its regime 
of oscillatory instability (cf. Fig. [jjj, for e — 0.5 and q = 5. We observe the configuration departing the steady state 
and giving rise to a persistent breathing state involving both components. 



VI. CONCLUSIONS 



In this paper, we have examined in detail the principal modes ari sing in the discrete vector lattices of AlGaAs 
waveguides, motivated by the recent experimental investigations of |23l l24j . We have constructed the full two- 
parameter bifurcation diagrams of all four principal branches of solutions. These include the transverse electric and 
transverse magnetic modes (occupying respectively each of the two components of the vector lattice), but also the 
linearly polarized (two components in or out of phase) and elliptically polarized modes (two components out of phase 
by odd multiples of 7r/2). The latter two modes, emerge from the former two in pairs through a pitchfork bifurcation 
that has been completely elucidated through the study of the system in the anti-continuum limit. In fact, the anti- 
continuum limit has provided with a very effective tool to analytically highlight existence and stability properties 
of the different branches as well as potential instabilities. We have continued this monoparametric picture (as a 
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FIG. 5: EP branch: The panels are analogous to those of Fig. [3] but for the EP branch; the top right panel, however, shows 
the non-vanishing power of the second component. Top left panels: q — 3.5. The middle and bottom panels correspond to 
e = 0.5 < €tm (left panels) and to e = 1 > e^M (right panels). In the latter case, the solution has already degenerated into 
a TM mode. Top right panels: q=4.5. The middle and bottom panels correspond to the stable case of e = 0.1 (left panels) 
and the unstable case of e = 1 (right panels). The bottom panel shows the two-parameter diagram for this mode. The region 
for q > 3 and < e < £tm> where this modes exists is separated by the dashed line into a stable and an unstable regime. For 
comparison the extension of the regions of stability/instability of the TM mode from Fig. |3]are also included. We note that the 
choice of star symbols (instead of circles as e.g. in Figs. |5|2J in the spatial profile of a n is used to indicate that the imaginary 
part of a n is shown. 




FIG. 6: This demonstrates the instability in the TE mode when e = 0.5 and q = 5.5. The left panel shows the space-time 
contour plot of the modulus square of the two fields (|a n (i)| 2 in the top panel and |&n(£)| 2 in the bottom panel). The right 
panel shows the evolution of the (square modulus of the) central site of the configuration as a function of time. 
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FIG. 7: An instability of the TM mode is demonstrated here with e = 0.5 and q — 3.5. The left panel again shows the 
space-time contour plot of the modulus square of the two fields (|a n (i)| 2 in the top panel and |6 n | 2 in the bottom panel). The 
right panel shows the evolution of the (square modulus of the) central site of the configuration as a function of time. 
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FIG. 8: The instability of a TM mode is shown here for e = 0.5 and q = 2.5, using exactly the same diagnostics as in Fig. Q 



function of the propagation constant of the solutions, but at zero coupling) to finite coupling fleshing out the regimes 
of stability and instability, and of existence as well as of disappearance of the relevant branches of solutions. The 
relevant two-parameter diagrams (Fig. |2]for the TE mode, Fig. |3|for the TM, Fig. Q]for the LP and Fig. [EJfor the EP 
mode) constitute the main result of the present paper. We have, however, complemented our existence and stability 
analysis with direct numerical "experiments" , illustrating the different instability scenario for the modes of interest 
and how they can either lead to the excitation of other branches, or even to the complete annihilation of the relevant 
localized waves. 





FIG. 9: Same as the previous 2 figures, but for the case of e = 1 and q = 2.5. Here the nonlinear mode is completely destroyed 
by the instability and results into small amplitude extended waves. 
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FIG. 10: An unstable elliptically polarized mode with e = 0.5 and q = 5 is monitored by the same diagnostics as in the previous 
figures. 

It is noteworthy that the process of obtaining such a multiparametric bifurcation picture has allowed us to reveal a 
considerably deeper understanding of the s yst em's parameter space, in comparison to earlier publications (such as e.g. 
the theoretical conclusions of the works of |23J,|2J]). ^ n particular, for instance, the limited parameter space examined 
in the latter setting (where a monoparametric continuation was used for e ~ 0.921 and apparently varying q only for 
q > 3) led the authors to infer that the TM mode is always unstable and similarly that the bifurcating EP mode is 
also unstable. Here, we have illustrated explicitly regions of stability of each of these modes in the parameter space 
of the system (which definitely exist for lower values of the coupling) . 

Finally, while the recent investigations of the various principal branches have offered a rather substantial under- 
standing of their properties, there are multiple directions of future interest in this topic that also appear within 
experimental reach. Such examples consist e.g. of the study of two-dimensional waveguide arrays and of discrete 
solitons [U and vortices |32l in them. Another such example is the study of more elaborate, multi-site constructions, 
such as twisted modes |33| even in one spatial dimension. While the above directions of examination of stationary 
states are of interest in their own right, an additional dimension of relevance revealed by our direct instability simula- 
tions of section V pertains to the systematic study of the existence and stability of breathing excitations due to their 
principal role in the system's dynamics. Such investigations are currently in progress and will be reported in future 
publications. 
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